clear; clc;
% 绘制膜电位 x(t) 与能量 H(t) 随时间的变化
% (a) f=0.154  的 x(t) 与 H(t)
% (b) f=0.3    的 x(t) 与 H(t)
% (c) f=0.7    的 x(t) 与 H(t)
%
% 其中 H(t) 根据论文给出的能量函数计算，

a     = 0.7; 
b     = 0.8; 
c     = 0.1; 
alpha = 0.25; 
beta  = 0.01; 
delta = 0.1;

% 外部激励幅值
A = 0.35;
f_list = [0.154, 0.3, 0.7];

y0 = [0.2; 0; 0.01];  % 初始条件

T_total = 1200;   % 总模拟时间
T_trans = 800;    % 舍去瞬态时间（后期数据用于绘图和计算平均能量）
h       = 0.01;   % RK4 积分步长

Purple = [0.5, 0, 0.5];   % 用于绘制膜电位 x(t)
Gold   = [1.0, 0.84, 0];    % 用于绘制能量 H(t)

% 3. 预分配存储数据
numCases = length(f_list);
t_all = cell(numCases,1);
x_all = cell(numCases,1);
H_all = cell(numCases,1);
meanH = zeros(numCases,1);

% 对每个频率进行积分并计算能量 H(t)
for i = 1:numCases
    f_val = f_list(i);
    params = struct('a', a, 'b', b, 'c', c, 'alpha', alpha, ...
                    'beta', beta, 'delta', delta, 'A', A, 'f', f_val);
                
    % 数值积分：使用自写 RK4 积分器
    [t, Y] = rk4(@(t,y) hybridNeuron(t, y, params), [0, T_total], y0, h);
    
    % 提取 x, y, z 数据
    x_data = Y(:,1);
    y_data = Y(:,2);
    z_data = Y(:,3);
    
    % 根据论文公式 (7) 计算能量 H(t)
    % H = 0.5*x^2 + (1/(2*c))*y^2 + 0.5*b*z*sin(beta*z)*y
    H_data = 0.5 * x_data.^2 + (1/(2*c)) * (y_data.^2) + 0.5 * b .* z_data .* sin(beta .* z_data) .* y_data;
    
    % 舍去瞬态数据
    idx = (t > T_trans);
    t_plot = t(idx);
    x_plot = x_data(idx);
    H_plot = H_data(idx);
    
    % 计算平均能量 <H>（取后期稳态数据均值）
    meanH(i) = mean(H_plot);
    
    % 保存数据
    t_all{i} = t_plot;
    x_all{i} = x_plot;
    H_all{i} = H_plot;
end

tlo = tiledlayout(3,2);

for i = 1:numCases
    % 左侧子图：绘制膜电位 x(t)
    nexttile(2*(i-1)+1);
    plot(t_all{i}, x_all{i}, 'Color', Purple, 'LineWidth', 1.5);
    xlabel('\tau'); ylabel('x');
    title(sprintf('f = %.3f', f_list(i)));
    grid on;
    
    % 右侧子图：绘制能量 H(t)
    nexttile(2*i);
    plot(t_all{i}, H_all{i}, 'Color', Gold, 'LineWidth', 1.5);
    xlabel('\tau'); ylabel('H');
    title(sprintf('f = %.3f; <H> = %.4f', f_list(i), meanH(i)));
    grid on;
end
